% $Id: adjoint.tex 204 2009-04-20 01:59:36Z jsibert $
%
% Author: David Fournier
% Copyright (c) 2008 Regents of the University of California
%

\mysection{The necessity for adjoint code}
When you write code for variable objects in \ADM\ all the
derivatives are calculated for you in a tranparent manner. 
to accomplish this \ADM\ must save certain information for
later use. We shall refer to this as derivative information.
Each arithmetic operation generates about 32 bytes of
derivative information.  If you have some simple function which
has 20 arithmetic operations it will therefore
generate 640 bytes of derivative information every time
it is called. The purpose of writing adjoint code is to
reduce the amount of derivative information which must
be calculated. For a function which is called many times
this can greatly reduce the amount of derivative information which
must be stored.

In this chapter we investiage how to write and debug
adjoint code. To begin we investigate how to write adjoint code for a
simple functioin which takes 1--4 independent variables and returns
one dependent variable. The adjoint code for such functions is
simpler to write than that for a general function such as the
singular value decomposition of a matrix which we will consider later.

\mysection{Writing Adjoint Code -- a simple case}

Consider a simple function $f$ which takes 1 independent variable $x$
and returns a dependent variable $y$. i.e.
$$y=f(x)$$
where $f(x)=\exp(-x^2/2)$. The code for this example can be written like
\beginexample
  dvariable errf(const prevariable& x)
  {
    return exp(-0.5*square(x));
  }
\endexample
\noindent
There are three arithmetic operations here, square, multiplication, and exponentiation
so that 96 bytes of derivative information will be generated. (Actually the return operation
also generates 32 bytes of derivative information but we will ignore that for now.)
A less efficient way to write the code (but more useful for 
showing adjoint code for this simple example) would be
\beginexample
  dvariable errf(const prevariable& x)
  {
    dvariable y;
    y=exp(-0.5*square(x));
    return y;
  }
\endexample
\noindent
Here is the same code with the derivative calculated by the one line of adjoint
code {\tt double dfx=-value(x)*value(y);}.
\beginexample
  dvariable errf(const prevariable& x)
  {
    dvariable y;
    value(y)=exp(-0.5*square(value(x)));
    double dfx=-value(x)*value(y);
    AD_SET_DERIVATIVES(y,x,dfx);  // 1 dependent variable 
    return y;
  }
\endexample
So what is going on here? Consider the line
    {\tt value(y)=exp(-0.5*square(value(x))); }
The {\tt value} function returns a constant type, that is, a double
which has the same value as the corresponding dvariable or prevariable.
In fact it is the {\bf same} object that is shares the same address,
but the type has been changed to double so that the above line of
code assigns the value {\tt exp(-0.5*square(value(x))); } to y
but without generating any derivative code. Simliarily sincd
the calculations are made on {\tt value(x)} these calculations
will not generate any derivative code. So it is the responsibility of the
programmer to calculate the derivative code and store it where it
can be used later. The line
    {\tt double dfx=-value(x)*value(y);} calculates the derivative $f^\prime(x)$
of $y$ with respect to $x$ and stores it with the line of code
    {\tt AD\_SET\_DERIVATIVES(y,x,dfx);} 
This code will only generate 32 bytes of derivative information.


\mysection{Debugging Adjoint Code -- a simple case}

The simplest way to debug the adjoint code is to put your new
function into an \ADM\ template file and use the {\tt -dd 1}
command line option to call the derivative checker.
\XX{-dd}{command line argument}
\XX{command line arguments}{-dd}
\beginexample
DATA_SECTION
PARAMETER_SECTION
  init_number x
 !! x=2;
  objective_function_value f;
PROCEDURE_SECTION
  f=square(errf(x);
  
GLOBALS_SECTION
  #include <admodel.h>
  dvariable errf(const prevariable& x)
  {
    dvariable y;
    value(y)=exp(-0.5*square(value(x)));
    double dfx=-value(x)*value(y);
    AD_SET_DERIVATIVES(y,x,dfx);  // 1 dependent variable 
    return y;
  }
\endexample
\mysection{Adjoint code for more than one indepdendent variable}

The following code shows how to write the adjoint code for a function with two
\beginexample
independent variables.
DATA_SECTION
  vector lengths(1,10)
  vector ages(1,10)
  !! lengths.fill_seqadd(1,1);
  !! ages.fill_seqadd(1,1);
  !! lengths=sqrt(lengths);

PARAMETER_SECTION
  init_bounded_number linf(0,10) 
  init_bounded_number rho(0,1) 
  objective_function_value f;
PROCEDURE_SECTION
  for (int i=1;i<=10;i++)
    f+=square(lengths(i)-vb_growth(linf,rho,ages(i)));
  
GLOBALS_SECTION
  #include <admodel.h>

  dvariable vb_growth(const prevariable& linf, const prevariable& rho, 
    double t)
  {
    double clinf=value(linf);
    double crho=value(rho);
    dvariable len;
    value(len)=clinf*(1-pow(crho,t));
    double dflinf=1-pow(crho,t);
    double dfrho=-clinf*t*pow(crho,t-1);
    AD_SET_DERIVATIVES2(len,rho,dfrho,linf,dflinf);  // 3 dependent variable 
    return len;
  }
\endexample
This approach to writing adjoint code has been implemented for function of
up to 4 independent variables.
\beginexample

DATA_SECTION
  vector lengths(1,10)
  vector ages(1,10)
  !! lengths.fill_seqadd(1,1);
  !! ages.fill_seqadd(1,1);
  !! lengths=sqrt(lengths);

PARAMETER_SECTION
  init_bounded_number linf(0,10) 
  init_bounded_number rho(0,1) 
  init_number t0 
  init_bounded_number gamma(.1,1.9) 
  objective_function_value f;
PROCEDURE_SECTION
  for (int i=1;i<=10;i++)
    f+=square(lengths(i)-vb_growth(linf,rho,t0,gamma,ages(i)));
  
GLOBALS_SECTION
  #include <admodel.h>

  dvariable vb_growth(const prevariable& linf, const prevariable& rho, 
    const prevariable& t0, const prevariable gamma,double t)
  {
    double clinf=value(linf);
    double ct0=value(t0);
    double crho=value(rho);
    double cgamma=value(gamma);
    dvariable len;
    value(len)=pow(clinf*(1-pow(crho,t-ct0)),cgamma);
    double tmp=cgamma*pow(clinf*(1-pow(crho,t-ct0)),cgamma-1);
    double dflinf=tmp*(1-pow(crho,t-ct0));
    double dft0=tmp*(clinf*log(crho)*pow(crho,t-ct0));
    double dfrho=-tmp*clinf*(t-ct0)*pow(crho,t-ct0-1);
    double dfgamma=value(len)*log(clinf*(1-pow(crho,t-ct0)));
    AD_SET_DERIVATIVES4(len,t0,dft0,rho,dfrho,linf,dflinf,gamma,dfgamma);  // 4 dependent variable 
    return len;
  }
\endexample

\mysection{Structured calculation of derivatives in adjoint code}
Until now we have deliberately calculated the derivatives with respect to the independent
variables in an {\it ad hoc} fashion. While this approach worsk for simple functions
it radidly becomes untenable when the functiuon is more complicated. In the following example
we have calculated the deivatives in a more structured fashion. Notice that to calculate the
derivatives every line of code in the function is repeated in the oposite
order (commented out of course) and the corresponding derivatves are calculated.

\beginexample

DATA_SECTION
  vector lengths(1,10)
  vector ages(1,10)
  !! lengths.fill_seqadd(1,1);
  !! ages.fill_seqadd(1,1);
  !! lengths=sqrt(lengths);

PARAMETER_SECTION
  init_bounded_number linf(0,10) 
  init_bounded_number rho(0,1) 
  init_number t0 
  init_bounded_number gamma(.1,1.9) 
  objective_function_value f;
PROCEDURE_SECTION
  for (int i=1;i<=10;i++)
    f+=square(lengths(i)-vb_growth(linf,rho,t0,gamma,ages(i)));
  
GLOBALS_SECTION
  #include <admodel.h>

  dvariable vb_growth(const prevariable& linf, const prevariable& rho, 
    const prevariable& t0, const prevariable gamma,double t)
  {
    double clinf=value(linf);
    double ct0=value(t0);
    double crho=value(rho);
    double cgamma=value(gamma);
    dvariable len;
    double u1=pow(crho,t-ct0);
    double u2=clinf*(1-u1); 
    value(len)=pow(u2,cgamma);
    double dflen=1.0;
    //value(len)=pow(u2,cgamma);
    double dfu2=dflen*cgamma*pow(u2,cgamma-1.0);
    double dfgamma=dflen*value(len)*log(u2);
    //double u2=clinf*(1-u1); 
    double dflinf=dfu2*(1-u1);
    double dfu1=-dfu2*clinf;
    //double u1=pow(crho,t-ct0);
    double dfrho=dfu1*(t-ct0)*pow(crho,t-ct0-1.0);
    double dft0=-dfu1*u1*log(crho);

    AD_SET_DERIVATIVES4(len,t0,dft0,rho,dfrho,linf,dflinf,gamma,dfgamma);  // 4 dependent variable 
    return len;
  }
\endexample
\mysection{General Adjoint Code}
So far the adjoint code has been for a simple function which has from 1 to 4
independent variables and returns one dependent variable. 
Now we condsider the general case where the function can take
any number of dependent variables and return any number of 
dependent variables and these variables can be in the form of numbers, 
vectors, or matrices.

\beginexample
 dvar_vector  operator * (_CONST dvar_matrix& m,_CONST dvar_vector& x )
 {
   if (x.indexmin() != m.colmin() || x.indexmax() != m.colmax())
   {
     cerr << " Incompatible array bounds in dvar_vector  operator * (_CONST dvar_matrix& m,_CONST dvar_vector& x)\n";
     ad_exit(21);
   }

   dvar_vector tmp(m.rowmin(),m.rowmax());
   double sum;

   for (int i=m.rowmin(); i<=m.rowmax(); i++)
   {
     sum=0.0;
     for (int j=x.indexmin(); j<=x.indexmax(); j++)
     {
       sum+=(m.elem(i)).elem_value(j)*x.elem_value(j);
     }
     tmp.elem_value(i)=sum;
   }
  save_identifier_string("PLACE4");
  x.save_dvar_vector_value();
  x.save_dvar_vector_position();
  save_identifier_string("PLACE3");
  m.save_dvar_matrix_value();
  m.save_dvar_matrix_position();
  save_identifier_string("PLACE2");
  tmp.save_dvar_vector_position();
  save_identifier_string("PLACE1");
  
  ADJOINT_CODE(dmdv_prod);
   return(tmp);
 }
\endexample
To calculate the adjoint code it will be necessary to have the
values of the matrix {\tt m} and the vector {\tt x}.
This is accompolshed with the
\beginexample
  x.save_dvar_vector_value();
  m.save_dvar_matrix_value();
\endexample
\noindent instructions.
also to calculate the derivatves it will be necessary to know where the
derivatives with repsect to the independent and dependent variables are 
located. This information is saved with the instructions
\beginexample
  x.save_dvar_vector_position();
  m.save_dvar_matrix_position();
  tmp.save_dvar_vector_position();
\endexample
\noindent Finally we need to save the name of the routine which calculates
the adjoint code so that it can be called at the appropriate time. To write the
code for the adjoint calculations keep in mind tha everything must be recovered
from the stack in the reverse order to which it was put on the stack.
this p;rocess can be a bit confusing and if you don't do it properly the 
stack will become corrupted and nothng will work. to help diagnose
problems function {\tt save\_identifier\_string} can be
used to put a string on the stack. This string value can be checked in the
adjoint code with the {\tt verify\_identifier\_string} function. At least two of
these functions should be left in any adjoint code so that stack
integrity can be monitored if problems show up later. For optimized code
they are not used and so contribute almost nothing to the overhead.

The adjoint code begins by reading the information which was saved on the stack.
An object of type {\tt dvar\_vector\_position} contains both the size and
address information associated with a {\tt dvar\_vector} which are needed to
recover or store derivative values or to build a {\tt dvector} with the
same shape as the {\t dvar\_vector}.
The function {\tt restore\_dvar\_vector\_derivatives} gets the values of the
derivatves with respect to the dependent variables sop that they can be used in the
adjoint code. the functions
\beginexample
   dfx.save_dvector_derivatives(x_pos);
   dfm.save_dmatrix_derivatives(m_pos);
\endexample
\noindent use the position information to save the derivatives with respect to the
independent variables in the appropriate places.
\beginexample
void dmdv_prod(void)
{
  verify_identifier_string("PLACE1");
  dvar_vector_position tmp_pos=restore_dvar_vector_position();
  verify_identifier_string("PLACE2");
  dvar_matrix_position m_pos=restore_dvar_matrix_position();
  dmatrix m=restore_dvar_matrix_value(m_pos);
  verify_identifier_string("PLACE3");
  dvar_vector_position x_pos=restore_dvar_vector_position();
  dvector x=restore_dvar_vector_value(x_pos);
  verify_identifier_string("PLACE4");
  dvector dftmp=restore_dvar_vector_derivatives(tmp_pos);
  
  dmatrix dfm(m_pos);
  dvector dfx(x_pos.indexmin(),x_pos.indexmax());
  dfm.initialize();
  dfx.initialize();
  
   double dfsum;
   for (int i=m.rowmax(); i>=m.rowmin(); i--)
   {
     // tmp.elem_value(i)=sum;
     dfsum=dftmp.elem(i);
     for (int j=x.indexmax(); j>=x.indexmin(); j--)
     {
       //sum+=(m.elem(i)).elem_value(j)*x.elem_value(j);
       dfm.elem(i,j)+=dfsum*x.elem(j);
       dfx.elem(j)+=dfsum*m.elem(i,j);
     }
     //sum=0.0;
     dfsum=0.0;
   }
   dfx.save_dvector_derivatives(x_pos);
   dfm.save_dmatrix_derivatives(m_pos);
}
\endexample
